#ifndef __CLASS_MATRIX_CPP
#define __CLASS_MATRIX_CPP
#include "class_ma_Complex.h"
//#include "class_maf.cpp"

// class_ma.cpp
/*inline float conj(float &a)
{
	return a;
}
double norm_inf(matrix_ &a)
{
	int i,j,I,J;
   double max=0,sum;
   I=a.rows;J=a.columns;
   for(i=0;i<I;i++)
   {
      sum=0;
   	for(j=0;j<J;j++)
      {
      	sum+=abs(a[i*J+j]);
      }
      if(sum>max)max=sum;
   }
   return max;
}
double norm_1(matrix_ &a)
{
	int i,j,I,J;
   double max=0,sum;
   I=a.rows;J=a.columns;
   for(j=0;j<J;j++)
   {
      sum=0;
   	for(i=0;i<I;i++)
      {
      	sum+=abs(a[i*J+j]);
      }
      if(sum>max)max=sum;
   }
   return max;
}

void CG1(float Er,int order,COMPLEX *A,COMPLEX *x,COMPLEX *b);
void matrix_::conjugate()
{
	int i,I;
	I=rows*columns;
	for(i=0;i<I;i++)
	array[i]=conj(array[i]);
};
void conj(matrix_ &a)
{
	int i,I;
	I=a.rows*a.columns;
	for(i=0;i<I;i++)
	a.array[i]=conj(a.array[i]);
} */
matrix_ matrix_::operator=(const COMPLEX &a)
{
	int i,I;
	I=rows*columns;
	for(i=0;i<I;i++)
	array[i]=a;
	return *this;
}

matrix_::matrix_(int row,int column)
{
  	int i,I;
   I=row*column;
   rows=row;
   columns=column;
   array=new Complex[I];
   if(array==NULL)
   {
   	cout<<"\nCan't allocate memory for the current matrix.";
      cin>>i;
   	exit(0);
   };
   for(i=0;i<I;i++)array[i]=COMPLEX(0);
   
}
/*
void change_rows_cols(matrix_ &A,int row,int column)
{
  	int i,I;
   I=row*column;
   A.rows=row;
   A.columns=column;

   for(i=0;i<I;i++)A.array[i]=COMPLEX(0);

   //sysout("\n construction");
}

void change_rows_cols_1(matrix_ &A,int row,int column)
{
	A.rows=row;
   A.columns=column;
   
   //sysout("\n construction");
}

matrix_ matrix_::slice(int row1,int col1,int row2,int col2)
{
	int col=columns;
   matrix_ m(row2,col2);
   COMPLEX *array=m.array;
   COMPLEX *array_this=array;
   for(int i=row1;i<row1+row2;i++)
   for(int j=col1;j<col1+col2;j++)
   {
   	array[(i-row1)*col2+(j-col1)]=array_this[i*col+j];
   };
   return(m);
}
matrix_ slice(const matrix_ &m1,int row1,int col1,int row2,int col2)
{
	int col=m1.columns;
   matrix_ m(row2,col2);
   COMPLEX *array=m.array;
   COMPLEX *array1=m1.array;
   for(int i=row1;i<row1+row2;i++)
   for(int j=col1;j<col1+col2;j++)
   {
   	array[(i-row1)*col2+(j-col1)]=array1[i*col+j];
   };
   return(m);
}

void slice(const matrix_ &m,const matrix_ &m1,int row1,int col1,int row2,int col2)
{
	int col=m1.columns;
   COMPLEX *array=m.array;
   COMPLEX *array1=m1.array;
   for(int i=row1;i<row1+row2;i++)
   for(int j=col1;j<col1+col2;j++)
   {
   	array[(i-row1)*col2+(j-col1)]=array1[i*col+j];
   };  
}

void set_value(const matrix_ &m,const matrix_ &m1,int row1,int col1,int row2,int col2)
{
	int col=m1.columns;
   COMPLEX *array=m.array;
   COMPLEX *array1=m1.array;
   for(int i=row1;i<row1+row2;i++)
   for(int j=col1;j<col1+col2;j++)
   {
   	array1[i*col+j]=array[(i-row1)*col2+(j-col1)];
   };  
}

matrix_ operator+(const matrix_ &m1,const matrix_ &m2)
{
      int i,r,c;
      r=m1.rows;
      c=m1.columns;
      matrix_ m(r,c);
   	//do the operations
      COMPLEX *array=m.array;
      COMPLEX *array1=m1.array;
      COMPLEX *array2=m2.array;
      for(i=0;i<(r*c);i++,array++,array1++,array2++)
      *array=*array1+*array2;
      return(m);
}
matrix_ operator-(const matrix_ &m1,const matrix_ &m2)
{
      int	 i,r,c;
      r=m1.rows;
      c=m1.columns;
      matrix_ m(r,c);
   	//do the operations
      COMPLEX *array=m.array;
      COMPLEX *array1=m1.array;
      COMPLEX *array2=m2.array;
      for(i=0;i<(r*c);i++,array++,array1++,array2++)
      *array=*array1-*array2;
      return(m);
}
matrix_ operator*(const matrix_ &m1,const matrix_ &m2)
{
      int	 i,j,k,tr,tc,c1;
      tr=m1.rows;
      tc=m1.columns;
      c1=m2.columns;
      matrix_ m(tr,c1);
      //do the operations
      COMPLEX *array=m.array;
      COMPLEX *array1=m1.array;
      COMPLEX *array2=m2.array;
      COMPLEX a,b;
      for(i=0;i<tr;i++)
      for(j=0;j<c1;j++)
      {
      	array[i*c1+j]=COMPLEX(0);
      	for(k=0;k<tc;k++)
      	{
      		array[i*c1+j]+=array1[i*tc+k]*array2[k*c1+j];

         };
      };
      return(m);
}
matrix_ operator*(const matrix_ &m1,COMPLEX x)
{
      int i,r,c;
      r=m1.rows;
      c=m1.columns;
      matrix_ m(r,c);
   	//do the operations
      COMPLEX *array=m.array;
      COMPLEX *array1=m1.array;
      for(i=0;i<(r*c);i++,array++,array1++)
      *array=*array1*x;
      return(m);
}
matrix_ operator*(COMPLEX x,const matrix_ &m1)
{
      int	 i,r,c;
      r=m1.rows;
      c=m1.columns;
      matrix_ m(r,c);
   	//do the operations
      COMPLEX *array=m.array;
      COMPLEX *array1=m1.array;
      for(i=0;i<(r*c);i++,array++,array1++)
      *array=*array1*x;
      return(m);
}
matrix_ operator/(const matrix_ &m1,COMPLEX x)
{
      int i,r,c;
      if(x==COMPLEX(0,0))
      {
      	cout<<"\nFAULT !!!Divid by zero.";
         cout<<"\nPress any key to end this game.";
      	int k;
         cin>>k;
         exit(0);};
      r=m1.rows;
      c=m1.columns;
      matrix_ m(r,c);
   	//do the operations
      COMPLEX *array=m.array;
      COMPLEX *array1=m1.array;
      for(i=0;i<(r*c);i++,array++,array1++)
      *array=*array1/x;
      return(m);
}

void invert(const matrix_ &m,int *pp,int *qq)
{
   float Error=1.e-24;
   COMPLEX z;
   int a,b;
   int i,j,k,l,mm,n;
	int ik,jk;
   float max;
   float ww;
   int mr=m.rows;
   int mc=m.columns;
   if(mr!=mc)
   {
   	cout<<"\nIt is not a n*n matrix";
      cout<<"the game is end!!input a number:";
      cin>>k;
      exit(0);
   };
   COMPLEX *array=m.array;

   COMPLEX temp;
   for(k=0;k<mr;k++)
   {
   	pp[k]=k;
      qq[k]=k;
   };
	for(k=0;k<mr;k++)
   {

   	ik=k;jk=k;
     	max=0.;
   	for(i=k;i<mr;i++)
      {
         mm=i*mc;
   		for(j=k;j<mc;j++)
   		{
				ww=abs(*(array+mm+j));
 	   		if(ww>max)
      		{
      			ik=i;
         		jk=j;
         		max=ww;
         	};
         }
      };
      if(max<=Error)
   	{
   		cout<<"\nFailue of inversion!!!";
	      cout<<"\nEncounted zero pricinple element.!!!";
   	   cout<<"\nInput a number, and then the game will be end!!!";
      	cout<<"\nInput a number:";cin>>k;
	      exit(0);
   	};
	   if(ik!=k)
   	{
         mm=k*mc;n=ik*mc;
   		for(j=0;j<mc;j++)
	      {
   	   	a=mm+j;
      	   b=n+j;
      		temp=*(array+a);
	         *(array+a)=*(array+b);
   	      *(array+b)=temp;
	      };
   	};
	   if(jk!=k)
   	{
         for(i=0;i<mr;i++)
	      {
   	   	a=i*mc+k;
      	   b=i*mc+jk;
         	temp=*(array+a);
	         *(array+a)=*(array+b);
   	      *(array+b)=temp;
	      };
   	};
	   *(pp+k)=ik;*(qq+k)=jk;
      mm=k*mc;
   	z=float(1.)/array[mm+k];
      for(j=0;j<mc;j++)array[mm+j]*=z;
      array[mm+k]=z;

      for(i=0;i<mr;i++)
      {
         mm=i*mc;n=k*mc;
      	if(i!=k)
         {
         	for(j=0;j<mc;j++)
	      	if(j!=k)*(array+mm+j)-=*(array+mm+k)**(array+n+j);
         	*(array+mm+k)*=(-z);
         };
   	};
   };
   for(l=0;l<mr-1;l++)
   {
   	k=mr-2-l;
   	ik=pp[k];
   	if(ik!=k)
      for(i=0;i<mr;i++)
      {
      	a=i*mc+k;
         b=i*mc+ik;
      	temp=array[a];
         array[a]=array[b];
         array[b]=temp;
      };
      jk=qq[k];
      if(jk!=k)
      {
         mm=k*mc;n=jk*mc;
      	for(j=0;j<mc;j++)
      	{
      		a=mm+j;
         	b=n+j;
      		temp=*(array+a);
         	*(array+a)=*(array+b);
         	*(array+b)=temp;
      	};
      }
   };
}

void invert(double *m,int mr,int *pp,int *qq)
{
   float Error=1.e-24;
   double z;
   int a,b;
   int i,j,k,l,mm,n;
	int ik,jk;
   float max;
   float ww;
   int mc=mr;

   double *array=m;

   double temp;
   for(k=0;k<mr;k++)
   {
   	pp[k]=k;
      qq[k]=k;
   };
	for(k=0;k<mr;k++)
   {

   	ik=k;jk=k;
     	max=0.;
   	for(i=k;i<mr;i++)
      {
         mm=i*mc;
   		for(j=k;j<mc;j++)
   		{
				ww=fabs(*(array+mm+j));
 	   		if(ww>max)
      		{
      			ik=i;
         		jk=j;
         		max=ww;
         	};
         }
      };
      if(max<=Error)
   	{
   		cout<<"\nFailue of inversion!!!";
	      cout<<"\nEncounted zero pricinple element.!!!";
   	   cout<<"\nInput a number, and then the game will be end!!!";
      	cout<<"\nInput a number:";cin>>k;
	      exit(0);
   	};
	   if(ik!=k)
   	{
         mm=k*mc;n=ik*mc;
   		for(j=0;j<mc;j++)
	      {
   	   	a=mm+j;
      	   b=n+j;
      		temp=*(array+a);
	         *(array+a)=*(array+b);
   	      *(array+b)=temp;
	      };
   	};
	   if(jk!=k)
   	{
         for(i=0;i<mr;i++)
	      {
   	   	a=i*mc+k;
      	   b=i*mc+jk;
         	temp=*(array+a);
	         *(array+a)=*(array+b);
   	      *(array+b)=temp;
	      };
   	};
	   *(pp+k)=ik;*(qq+k)=jk;
      mm=k*mc;
   	z=float(1.)/array[mm+k];
      for(j=0;j<mc;j++)array[mm+j]*=z;
      array[mm+k]=z;

      for(i=0;i<mr;i++)
      {
         mm=i*mc;n=k*mc;
      	if(i!=k)
         {
         	for(j=0;j<mc;j++)
	      	if(j!=k)*(array+mm+j)-=*(array+mm+k)**(array+n+j);
         	*(array+mm+k)*=(-z);
         };
   	};
   };
   for(l=0;l<mr-1;l++)
   {
   	k=mr-2-l;
   	ik=pp[k];
   	if(ik!=k)
      for(i=0;i<mr;i++)
      {
      	a=i*mc+k;
         b=i*mc+ik;
      	temp=array[a];
         array[a]=array[b];
         array[b]=temp;
      };
      jk=qq[k];
      if(jk!=k)
      {
         mm=k*mc;n=jk*mc;
      	for(j=0;j<mc;j++)
      	{
      		a=mm+j;
         	b=n+j;
      		temp=*(array+a);
         	*(array+a)=*(array+b);
         	*(array+b)=temp;
      	};
      }
   };
}
void invert(double *m,int mr)
{
   float Error=1.e-24;
   double z;
   int a,b;
   int i,j,k,l,mm,n;
	int ik,jk;
   int *qq,*pp;
   float max;
   float ww;
   int mc=mr;

   double *array=m;

   double temp;

   pp=new int[mr];
   qq=new int[mc];
   for(k=0;k<mr;k++)
   {
   	pp[k]=k;
      qq[k]=k;
   };
	for(k=0;k<mr;k++)
   {

   	ik=k;jk=k;
     	max=0.;
   	for(i=k;i<mr;i++)
      {
         mm=i*mc;
   		for(j=k;j<mc;j++)
   		{
				ww=fabs(*(array+mm+j));
 	   		if(ww>max)
      		{
      			ik=i;
         		jk=j;
         		max=ww;
         	};
         }
      };
      if(max<=Error)
   	{
   		cout<<"\nFailue of inversion!!!";
	      cout<<"\nEncounted zero pricinple element.!!!";
   	   cout<<"\nInput a number, and then the game will be end!!!";
      	cout<<"\nInput a number:";cin>>k;
	      exit(0);
   	};
	   if(ik!=k)
   	{
         mm=k*mc;n=ik*mc;
   		for(j=0;j<mc;j++)
	      {
   	   	a=mm+j;
      	   b=n+j;
      		temp=*(array+a);
	         *(array+a)=*(array+b);
   	      *(array+b)=temp;
	      };
   	};
	   if(jk!=k)
   	{
         for(i=0;i<mr;i++)
	      {
   	   	a=i*mc+k;
      	   b=i*mc+jk;
         	temp=*(array+a);
	         *(array+a)=*(array+b);
   	      *(array+b)=temp;
	      };
   	};
	   *(pp+k)=ik;*(qq+k)=jk;
      mm=k*mc;
   	z=float(1.)/array[mm+k];
      for(j=0;j<mc;j++)array[mm+j]*=z;
      array[mm+k]=z;

      for(i=0;i<mr;i++)
      {
         mm=i*mc;n=k*mc;
      	if(i!=k)
         {
         	for(j=0;j<mc;j++)
	      	if(j!=k)*(array+mm+j)-=*(array+mm+k)**(array+n+j);
         	*(array+mm+k)*=(-z);
         };
   	};
   };
   for(l=0;l<mr-1;l++)
   {
   	k=mr-2-l;
   	ik=pp[k];
   	if(ik!=k)
      for(i=0;i<mr;i++)
      {
      	a=i*mc+k;
         b=i*mc+ik;
      	temp=array[a];
         array[a]=array[b];
         array[b]=temp;
      };
      jk=qq[k];
      if(jk!=k)
      {
         mm=k*mc;n=jk*mc;
      	for(j=0;j<mc;j++)
      	{
      		a=mm+j;
         	b=n+j;
      		temp=*(array+a);
         	*(array+a)=*(array+b);
         	*(array+b)=temp;
      	};
      }
   };
   delete(qq);
   delete(pp);
}*/
void invert(const matrix_ &m)
{
   double Error=1.e-24;
   Complex z;
   int a,b;
   int i,j,k,l,mm,n;
	int ik,jk;
   int *pp,*qq;
   double max;
   double ww;
   int mr=m.rows;
   int mc=m.columns;
   if(mr!=mc)
   {
   	cout<<"\nIt is not a n*n matrix";
      cout<<"the game is end!!input a number:";
      cin>>k;
      exit(0);
   };
   Complex *array=m.array;

   Complex temp;
   pp=new int[mr];
   qq=new int[mc];
   for(k=0;k<mr;k++)
   {
   	pp[k]=k;
      qq[k]=k;
   };
	for(k=0;k<mr;k++)
   {

   	ik=k;jk=k;
     	max=0.;
   	for(i=k;i<mr;i++)
      {
         mm=i*mc;
   		for(j=k;j<mc;j++)
   		{
				ww=abs(*(array+mm+j));
 	   		if(ww>max)
      		{
      			ik=i;
         		jk=j;
         		max=ww;
         	};
         }
      };
      if(max<=Error)
   	{
   		cout<<"\nFailue of inversion!!!";
	      cout<<"\nEncounted zero pricinple element.!!!";
   	   cout<<"\nInput a number, and then the game will be end!!!";
      	cout<<"\nInput a number:";cin>>k;
	      exit(0);
   	};
	   if(ik!=k)
   	{
         mm=k*mc;n=ik*mc;
   		for(j=0;j<mc;j++)
	      {
   	   	a=mm+j;
      	   b=n+j;
      		temp=*(array+a);
	         *(array+a)=*(array+b);
   	      *(array+b)=temp;
	      };
   	};
	   if(jk!=k)
   	{
         for(i=0;i<mr;i++)
	      {
   	   	a=i*mc+k;
      	   b=i*mc+jk;
         	temp=*(array+a);
	         *(array+a)=*(array+b);
   	      *(array+b)=temp;
	      };
   	};
	   *(pp+k)=ik;*(qq+k)=jk;
      mm=k*mc;
   	z=double(1.)/array[mm+k];
      for(j=0;j<mc;j++)array[mm+j]*=z;
      array[mm+k]=z;

      for(i=0;i<mr;i++)
      {
         mm=i*mc;n=k*mc;
      	if(i!=k)
         {
         	for(j=0;j<mc;j++)
	      	if(j!=k)*(array+mm+j)-=*(array+mm+k)**(array+n+j);
         	*(array+mm+k)*=(-z);
         };
   	};
   };
   for(l=0;l<mr-1;l++)
   {
   	k=mr-2-l;
   	ik=pp[k];
   	if(ik!=k)
      for(i=0;i<mr;i++)
      {
      	a=i*mc+k;
         b=i*mc+ik;
      	temp=array[a];
         array[a]=array[b];
         array[b]=temp;
      };
      jk=qq[k];
      if(jk!=k)
      {
         mm=k*mc;n=jk*mc;
      	for(j=0;j<mc;j++)
      	{
      		a=mm+j;
         	b=n+j;
      		temp=*(array+a);
         	*(array+a)=*(array+b);
         	*(array+b)=temp;
      	};
      }
   };
   delete [] qq;
   delete [] pp;
}
/*
matrix_ matrix_::invert(float Er)
{
    int i,j,k;
    int mr=rows;
    int mc=columns;
    matrix_ m(mr,mc);
    COMPLEX *b,*x;
    if(mr!=mc)
    {
    	cout<<"\nIt is not a n*n matrix";
      cout<<"the game is end!!input a number:";
      cin>>k;
      exit(0);
    };
    b=new COMPLEX[mr];
    x=new COMPLEX[mr];
    COMPLEX *array=array;
    COMPLEX *marray=m.array;
    for(j=0;j<mc;j++)
    {
    	for(i=0;i<mr;i++)
      b[i]=COMPLEX(0);
      b[j]=COMPLEX(1);
      CG1(Er,mr,array,x,b);
      //cout<<"end of CG1.";int key;cin>>key;
      for(i=0;i<mr;i++)
      marray[i*mc+j]=x[i];
    };
    delete b;
    delete x;
    return m;
}
void CG1(float Er,int order,COMPLEX *A,COMPLEX *x,COMPLEX *b)
{
	int i,j,k,step;
	COMPLEX *r,*r1,*p,*w,*v,*c,a,a2,Ru0,Ruk,Ruk1,Beta;
   	float rr,value;
	r=new COMPLEX[order];
   r1=new COMPLEX[order];
	p=new COMPLEX[order];
	w=new COMPLEX[order];
	v=new COMPLEX[order];
	c=new COMPLEX[order];
  	for(i=0;i<order;i++)x[i]=COMPLEX(1);
	//step 1.
	for(i=0;i<order;i++)
	{
		r[i]=b[i];
		for(j=0;j<order;j++)
		r[i]=r[i]-A[i*order+j]*x[j];
	};
	step=1;
	//cout<<"\nstep of iteration is "<<step;
	for(i=0;i<order;i++)p[i]=COMPLEX(0);
	for(j=0;j<order;j++)
	for(i=0;i<order;i++)
	p[i]=p[i]+A[j*order+i]*r[j];
	Ru0=COMPLEX(0);
	for(i=0;i<order;i++)Ru0=Ru0+p[i]*p[i];
	k=0;
	Ruk=Ru0;
   do{
   		do{
					k=k+1;
					//step 2.
					for(i=0;i<order;i++)w[i]=COMPLEX(0);
					for(i=0;i<order;i++)
					for(j=0;j<order;j++)w[i]=w[i]+A[i*order+j]*p[j];
					a2=COMPLEX(0);
					for(i=0;i<order;i++)a2=a2+w[i]*w[i];
					a=Ruk/a2;
               for(i=0;i<order;i++)x[i]=x[i]+a*p[i];
  					for(i=0;i<order;i++)
					for(j=0;j<order;j++)r[i]=r[i]-a*A[i*order+j]*p[j];
					for(i=0;i<order;i++)v[i]=COMPLEX(0);
               for(j=0;j<order;j++)
					for(i=0;i<order;i++)
					v[i]=v[i]+A[j*order+i]*r[j];
					Ruk1=COMPLEX(0);
					for(i=0;i<order;i++)Ruk1=Ruk1+v[i]*v[i];
					Beta=Ruk1/Ruk;
               for(i=0;i<order;i++)p[i]=v[i]+Beta*p[i];
					//step 3.
					Ruk=Ruk1;
					step+=1;
				  	//cout<<"\nstep of iteration is "<<step;
				   //cout<<"\nRuk1="<<Ruk1<<"Ru0*Er="<<Ru0*Er;
    			}while(abs(Ruk1)>=abs(Ru0*Er));
    		rr=0;
    		for(i=0;i<order;i++)
    		{
    			r1[i]=COMPLEX(0);
     			for(j=0;j<order;j++)
     			r1[i]+=A[i*order+j]*x[j];
      		r1[i]=b[i]-r1[i];
    		};
    		value=0;
         for(i=0;i<order;i++){rr+=abs(r1[i]);value+=abs(x[i]);};
	}while((rr/value)>=Er);

	delete(r);
  	delete(r1);
	delete(p);
	delete(w);
	delete(v);
	delete(c);
}*/
inline matrix_::operator Complex *()
{
	return array;
}
/*matrix_ matrix_::compose(matrix_ &m2,int select)
{
   int i,j,k;
	int r1=rows;
   int c1=columns;
   int r2=m2.rows;
   int c2=m2.columns;
   int r;
   int c;
   matrix_ m1=*this;
   matrix_ m;
   if(select==1)
   {
      if(r1!=r2)
      {
      	cout<<"\nThe two matrixs' rows are not equal to eachother";
         cout<<"The game is end.Input a number:";cin>>r;
         exit(0);
      };
   	r=r1;
      c=c1+c2;
      m=matrix_(r,c);
      for(i=0;i<r;i++)
      {
      	for(j=0;j<c1;j++)
         m[i*c+j]=m1[i*c1+j];
         for(j=0;j<c2;j++)
         m[i*c+c1+j]=m2[i*c2+j];
      };
   }
   else if(select==2)
   {
      if(c1!=c2)
      {
      	cout<<"\nThe two matrixs' cols are not equal to eachother";
         cout<<"The game is end.Input a number:";cin>>r;
         exit(0);
      };
   	r=r1+r2;
      c=c1;
      m=matrix_(r,c);
      for(i=0;i<r1;i++)
      for(j=0;j<c;j++)
      m[i*c+j]=m1[i*c+j];
      k=r1*c;
      for(i=0;i<r2;i++)
      for(j=0;j<c;j++)
      m[k+i*c+j]=m2[i*c+j];
   };
   return(m);
}

matrix_ compose(matrix_ &m1,matrix_ &m2,int select)
{
   int i,j,k;
	int r1=m1.rows;
   int c1=m1.columns;
   int r2=m2.rows;
   int c2=m2.columns;
   int r;
   int c;
   matrix_ m;
   if(select==1)
   {
      if(r1!=r2)
      {
      	cout<<"\nThe two matrixs' rows are not equal to eachother";
         cout<<"The game is end.Input a number:";cin>>r;
         exit(0);
      };
   	r=r1;
      c=c1+c2;
      m=matrix_(r,c);
      for(i=0;i<r;i++)
      {
      	for(j=0;j<c1;j++)
         m[i*c+j]=m1[i*c1+j];
         for(j=0;j<c2;j++)
         m[i*c+c1+j]=m2[i*c2+j];
      };
   }
   else if(select==2)
   {
      if(c1!=c2)
      {
      	cout<<"\nThe two matrixs' cols are not equal to eachother";
         cout<<"The game is end.Input a number:";cin>>r;
         exit(0);
      };
   	r=r1+r2;
      c=c1;
      m=matrix_(r,c);
      for(i=0;i<r1;i++)
      for(j=0;j<c;j++)
      m[i*c+j]=m1[i*c+j];
      k=r1*c;
      for(i=0;i<r2;i++)
      for(j=0;j<c;j++)
      m[k+i*c+j]=m2[i*c+j];
   };
   return(m);
}


//This is the modified compose function, All the array operations are
//replaced by the pointer opeations
void compose(matrix_ &mm,matrix_ &mm1,matrix_ &mm2,int select)
{
   int i,j;
   COMPLEX *m=mm.array,*m1=mm1.array,*m2=mm2.array;
	int r1=mm1.rows;
   int c1=mm1.columns;
   int r2=mm2.rows;
   int c2=mm2.columns;

   if(select==1)
   {
      if(r1!=r2)
      {
      	cout<<"\nThe two matrixs' rows are not equal to eachother";
         cout<<"The game is end.Input a number:";cin>>r1;
         exit(0);
      };
      for(i=0;i<r1;i++)
      {
      	for(j=0;j<c1;j++)
         {
         	*m=*m1;
            m++;
            m1++;
         }
         for(j=0;j<c2;j++)
         {
         	*m=*m2;
            m++;
            m2++;
         }
      };
   }
   else if(select==2)
   {
      if(c1!=c2)
      {
      	cout<<"\nThe two matrixs' cols are not equal to eachother";
         cout<<"The game is end.Input a number:";cin>>r1;
         exit(0);
      };

      for(i=0;i<r1;i++)
      for(j=0;j<c1;j++)
      {
      	*m=*m1;
         m++;
         m1++;
      }
      for(i=0;i<r2;i++)
      for(j=0;j<c2;j++)
      {
         *m=*m2;
         m++;
         m2++;
      }
   };
}



matrix_ matrix_::multi(const matrix_ &a)
{
	matrix_ m;
   cout<<"\nbegin multi";int k;cin>>k;
   m=*this*a;
   cout<<"\nend multi";cin>>k;
   return m;
}
void matrix_::output()
{
	int i,j;
	int mr,mc;
	mr=rows;
	mc=columns;
	for(i=0;i<mr;i++)
	{
		cout<<"\n";
		for(j=0;j<mc;j++)
		{
			cout.setf(ios::scientific);
			cout<<" "<<array[i*mc+j];
			if(j%6==0)cout<<"\n";
			cout.unsetf(ios::scientific);
		}
	};
}
matrix_ operator-(matrix_ &m)
{
   int i;
	int mr=m.rows;
   int mc=m.columns;
   int mrc=mr*mc;
   matrix_ m1(mr,mc);
   for(i=0;i<mrc;i++)
   m1[i]=-m[i];
   return(m1);
}
matrix_ operator+(matrix_ &m)
{
	return(m);
}
void conjugate(matrix_ &m)
{
	int i;
	int mr=m.rows;
   int mc=m.columns;
   int mrc=mr*mc;

   for(i=0;i<mrc;i++)
   m[i]=conj(m[i]);

}	

void foutput(matrix_ &a)
{
	ofstream fout;
	char name[20];
	cout<<"\nPlease input the file name you want to store your file to!";
	cin>>name;
	fout.open(name);
	int mr=a.rows;
	int mc=a.columns;
	int i,j;
	fout.setf(ios::scientific);
	for(i=0;i<mr;i++)
   {
		for(j=0;j<mc;j++)
		{
			fout<<real(a[i*mc+j])<<" ";
      }
      fout<<"\n";
	};
   for(i=0;i<mr;i++)
   {
		for(j=0;j<mc;j++)
		{
			fout<<imag(a[i*mc+j])<<" ";
      }
      fout<<"\n";
	};
	fout.unsetf(ios::scientific);
	fout.close();
}

void foutput1(matrix_ &a,char *name)
{
	ofstream fout;
	fout.open(name);
	int mr=a.rows;
	int mc=a.columns;
	int i,j;
	fout.setf(ios::scientific);
	for(i=0;i<mr;i++)
   {
		for(j=0;j<mc;j++)
		{
			fout<<real(a[i*mc+j])<<' ';
		};
      fout<<'\n';
   }
   for(i=0;i<mr;i++)
   {
		for(j=0;j<mc;j++)
		{
			fout<<imag(a[i*mc+j])<<' ';
		};
      fout<<'\n';
   }
   fout.unsetf(ios::scientific);
	fout.close();
}

void foutput(matrix_ &a,char *name)
{
	ofstream fout;
	fout.open(name);
	int mr=a.rows;
	int mc=a.columns;
	int i,j;
	fout.setf(ios::scientific);
	for(i=0;i<mr;i++)
	for(j=0;j<mc;j++)
	{
		fout<<real(a[i*mc+j])<<","<<imag(a[i*mc+j])<<"\n";
	};
	fout.unsetf(ios::scientific);
	fout.close();
}
void foutput(float *a,char *d,float *b,int n,char *name)
{
	ofstream fout;
	fout.open(name);
	int i;
	fout.setf(ios::scientific);
	for(i=0;i<n;i++)
	{
		fout<<a[i]<<d<<b[i]<<"\n";
	};
	fout.unsetf(ios::scientific);
	fout.close();
}

void foutputabs(matrix_ &a,char *name)
{
	ofstream fout;
	fout.open(name);
	int mr=a.rows;
	int mc=a.columns;
	int i,j;
	fout.setf(ios::scientific);
	for(i=0;i<mr;i++)
	{
		fout<<"\n";
		for(j=0;j<mc;j++)
		{
         if(abs(a[i*mc+j])>1e-6)
         fout<<1<<' ';
         else
         fout<<0<<' ';
			//fout<<abs(a[i*mc+j])<<' ';
		};
	};
	fout.unsetf(ios::scientific);
	fout.close();
}
void foutput_ma(matrix_ &a,char *name)
{
	ofstream fout;
	fout.open(name);
	int mr=a.rows;
	int mc=a.columns;
	int i,j;
	fout.setf(ios::scientific);
	for(i=0;i<mr;i++)
	{
		fout<<"\n";
		for(j=0;j<mc;j++)
		{
			fout<<a[i*mc+j]<<' ';

		};
	};
	fout.unsetf(ios::scientific);
	fout.close();
}
void unite_judge(matrix_ &a,char *name)
{
	ofstream fout;
	fout.open(name,ios::ate);
	int mr=a.rows;
	int mc=a.columns;
	int i,j;
	fout.setf(ios::scientific);

   fout<<"\n\n\nRows="<<a.rows;
   fout<<"\ncols="<<a.columns;

	for(i=0;i<mr;i++)
	{
		for(j=0;j<mc;j++)
		{
			if(i!=j&&abs(a[i*mc+j])>1e-4)
         {
         	fout<<"\nA["<<i<<"]["<<j<<"]="<<a[i*mc+j]<<"............."<<abs(a[i*mc+j]);
         }
         else if(i==j&&abs(a[i*mc+j]-1)>1e-4)
         {
            fout<<"\nA["<<i<<"]["<<j<<"]="<<a[i*mc+j]<<"............."<<abs(a[i*mc+j]);
         }

		};
	};
	fout.unsetf(ios::scientific);
	fout.close();
}
void finput(matrix_ &a,int mr,int mc)
{
	ifstream fin;
	char name[20],d;
	cout<<"\nPlease input the file name you want to get data from!";
	cin>>name;
	fin.open(name);
	int i,j;
	float re,im;
	fin.setf(ios::scientific);
	for(i=0;i<mr;i++)
	for(j=0;j<mc;j++)
	{
		fin>>re>>d>>im;
		a[i*mc+j]=COMPLEX(re,im);
	}
	fin.unsetf(ios::scientific);
	fin.close();
}
void finput(matrix_ &a,int mr,int mc,char *name)
{
	ifstream fin;
	char d;
	fin.open(name);
	int i,j;
	float re,im;
	fin.setf(ios::scientific);
	for(i=0;i<mr;i++)
	for(j=0;j<mc;j++)
	{
		fin>>re>>d>>im;
		a[i*mc+j]=COMPLEX(re,im);
	}
	fin.unsetf(ios::scientific);
	fin.close();
}

void finput(matrix_ &a,char *name)
{
	ifstream fin;
	char d;
	fin.open(name);
	int i,j;
	float re,im;
   int mr=a.rows,mc=a.columns;
	fin.setf(ios::scientific);
	for(i=0;i<mr;i++)
	for(j=0;j<mc;j++)
	{
		fin>>re>>d>>im;
		a[i*mc+j]=COMPLEX(re,im);
	}
	fin.unsetf(ios::scientific);
	fin.close();
}

matrix_ H(const matrix_ &m)
{
	int i,j;
	int row=m.rows;
	int col=m.columns;
	matrix_ mt(col,row);
	for(i=0;i<row;i++)
	for(j=0;j<col;j++)
	mt.array[j*row+i]=conj(m.array[i*col+j]);
	return mt;
}
matrix_ T(const matrix_ &m)
{
	int i,j;
	int row=m.rows;
	int col=m.columns;
	matrix_ mt(col,row);
	for(i=0;i<row;i++)
	for(j=0;j<col;j++)
	mt.array[j*row+i]=m.array[i*col+j];
	return mt;
}
float abs_matrix(matrix_ &m)
{
	COMPLEX a=0;
	int i,I=m.rows*m.columns;
	for(i=0;i<I;i++)
	a+=m[i]*conj(m[i]);
	return sqrt(abs(a));
}


void free_ma(matrix_ &a)
{
	if(a.array!=NULL)
	{
		delete []a.array;
		a.array=NULL;
	};
	//delete a.p;
	
};

void copy(matrix_ &a,matrix_ &b)
{
	int i,terms=a.rows*a.columns;
   COMPLEX *aa=a.array,*bb=b.array;
	for(i=0;i<terms;i++)
	{
		*aa=*bb;
      aa++;bb++;
	};
}

//For some case that a and b have different number of entries
void copy(matrix_ &a,matrix_ &b,int terms)
{
	int i;
	for(i=0;i<terms;i++)
	{
		a[i]=b[i];
	};
}
void copyHforvector(matrix_ &a,matrix_ &b)
{
	int i,terms=a.rows*a.columns;
	//here suppose all of the two matrixs are vectors.
	for(i=0;i<terms;i++)
	{
		a[i]=conj(b[i]);
	};
}

//This is the copy one:
void copy(matrix_ &a,COMPLEX b)
{
	int i,terms=a.rows*a.columns;
   COMPLEX *aa=a.array;
	for(i=0;i<terms;i++)
	{
		*aa=b;
      aa++;
	};
}


//These are further modified add, sub, multi, Hmulti, and multiH.
//No array is used. All are substituted by the form of pointers.
void add(matrix_ &c,matrix_ &a,matrix_ &b)
{
	int i,I=a.rows*a.columns;
   COMPLEX *aa=a.array,*bb=b.array,*cc=c.array;
	for(i=0;i<I;i++)
   {
		*(cc+i)=*(aa+i)+*(bb+i);
   }
};
void sub(matrix_ &c,matrix_ &a,matrix_ &b)
{
	int i,I=a.rows*a.columns;
   COMPLEX *aa=a.array,*bb=b.array,*cc=c.array;
	for(i=0;i<I;i++)
	*(cc+i)=*(aa+i)-*(bb+i);
};*/
void multi(matrix_ &c,matrix_ &a,matrix_ &b)
{
	int i,j,k,I=a.rows,J=b.columns,K=a.columns,na,nb,IJ=I*J;
   Complex *aa=a.array,*bb=b.array,*cc=c.array;

   for(i=0;i<IJ;i++)
   *(cc+i)=Complex(0);

   na=0;
   for(i=0;i<I;i++)
   {
   	nb=0;
      for(k=0;k<K;k++)
      {
         for(j=0;j<J;j++)
         {
            *(cc+j)+=*(aa+na)**(bb+nb);
         	nb++;
         }
         na++;
      }
      cc+=J;
   }
}
/*void multi(matrix_ &c,matrixf_ &a,matrixf_ &b)
{
	int i,j,k,I=a.rows,J=b.columns,K=a.columns,na,nb,IJ=I*J;
   float *aa=a.array,*bb=b.array;
   COMPLEX *cc=c.array;

   for(i=0;i<IJ;i++)
   *(cc+i)=COMPLEX(0);

   na=0;
   for(i=0;i<I;i++)
   {
   	nb=0;
      for(k=0;k<K;k++)
      {
         for(j=0;j<J;j++)
         {
            *(cc+j)+=COMPLEX(*(aa+na)**(bb+nb));
         	nb++;
         }
         na++;
      }
      cc+=J;
   }
}
void multi(matrix_ &c,matrix_ &a,matrixf_ &b)
{
	int i,j,k,I=a.rows,J=b.columns,K=a.columns,na,nb,IJ=I*J;
   COMPLEX *aa=a.array,*cc=c.array;
   FLOAT *bb=b.array;
   for(i=0;i<IJ;i++)
   *(cc+i)=COMPLEX(0);

   na=0;
   for(i=0;i<I;i++)
   {
   	nb=0;
      for(k=0;k<K;k++)
      {
         for(j=0;j<J;j++)
         {
            *(cc+j)+=*(aa+na)**(bb+nb);
         	nb++;
         }
         na++;
      }
      cc+=J;
   }
}
void multi(matrix_ &c,matrixf_ &a,matrix_ &b)
{
	int i,j,k,I=a.rows,J=b.columns,K=a.columns,na,nb,IJ=I*J;
   FLOAT *aa=a.array;
   COMPLEX *bb=b.array,*cc=c.array;

   for(i=0;i<IJ;i++)
   *(cc+i)=COMPLEX(0);

   na=0;
   for(i=0;i<I;i++)
   {
   	nb=0;
      for(k=0;k<K;k++)
      {
         for(j=0;j<J;j++)
         {
            *(cc+j)+=*(aa+na)**(bb+nb);
         	nb++;
         }
         na++;
      }
      cc+=J;
   }
}
void Tmulti(matrix_ &c,matrix_ &a,matrix_ &b)
{
	int i,j,k,I=a.columns,J=b.columns,K=a.rows,IJ=I*J,na,nc;
   COMPLEX *aa=a.array,*bb=b.array,*cc=c.array;

   for(i=0;i<IJ;i++)
   *(cc+i)=COMPLEX(0);

   na=0;
   for(k=0;k<K;k++)
   {
   	nc=0;
   	for(i=0;i<I;i++)
      {
      	for(j=0;j<J;j++)
         {
            *(cc+nc)+=*(aa+na)**(bb+j);
         	nc++;
         }
         na++;
      }
      bb+=J;
   }
}

//c=H(a)*b;
void Hmulti(matrix_ &c,matrix_ &a,matrix_ &b)
{
	int i,j,k,I=a.columns,J=b.columns,K=a.rows,IJ=I*J,na,nc;
   COMPLEX *aa=a.array,*bb=b.array,*cc=c.array;

   for(i=0;i<IJ;i++)
   *(cc+i)=COMPLEX(0);

   na=0;
   for(k=0;k<K;k++)
   {
   	nc=0;
   	for(i=0;i<I;i++)
      {
      	for(j=0;j<J;j++)
         {
            *(cc+nc)+=conj(*(aa+na))**(bb+j);
         	nc++;
         }
         na++;
      }
      bb+=J;
   }
}

//c=a*H(b);
void multiH(matrix_ &c,matrix_ &a,matrix_ &b)
{
   int i,j,k,I=a.rows,J=b.rows,K=a.columns,IJ=I*J,nb,nc;
   COMPLEX *aa=a.array,*bb=b.array,*cc=c.array;

   for(i=0;i<IJ;i++)
   *(cc+i)=COMPLEX(0);

   nc=0;
   for(i=0;i<I;i++)
   {
      nb=0;
      for(j=0;j<J;j++)
      {
      	for(k=0;k<K;k++)
         {
         	*(cc+nc)+=*(aa+k)*conj(*(bb+nb));
            nb++;
         }
         nc++;
      }
      aa+=K;
   }
}

//Modified:
COMPLEX square(matrix_ &AA)
{
	COMPLEX f(0,0),*A=AA.array;


	int i,I=AA.rows*AA.columns;
	for(i=0;i<I;i++)
   {
		f+=conj(*A)**A;
      A++;
   }
	return f;
}

void square(COMPLEX &f,matrix_ &A)
{
	f=COMPLEX(0);

	int i,I=A.rows*A.columns;
	for(i=0;i<I;i++)
	f+=conj(A[i])*A[i];

}

//Modified:
void divid(matrix_ &cc,matrix_ &bb,COMPLEX &a)
{
   COMPLEX *c=cc.array,*b=bb.array;
	int i,I=cc.rows*cc.columns;
	if(abs(a)<1e-30)
	{
		cout<<"error in matrix operation c=b/a, a=COMPLEX(0)";
		cin>>i;
		exit(0);
	};
	for(i=0;i<I;i++)
	*(c+i)=*(b+i)/a;
}

//Modified:
void divid(matrix_ &cc,matrix_ &bb,Complex &a)
{
   COMPLEX *c=cc.array,*b=bb.array;
	int i,I=cc.rows*cc.columns;
	if(abs(a)<1e-30)
	{
		cout<<"error in matrix operation c=b/a, a=COMPLEX(0)";
		cin>>i;
		exit(0);
	};
	for(i=0;i<I;i++)
	*(c+i)=COMPLEX(Complex(*(b+i))/a);
}

//Modified:
void divid(matrix_ &cc,matrix_ &bb,float &a)
{
   COMPLEX *c=cc.array,*b=bb.array;
	int i,I=cc.rows*cc.columns;
	if(fabs(a)<1e-30)
	{
		cout<<"error in matrix operation c=b/a, a=COMPLEX(0)";
		cin>>i;
		exit(0);
	};
	for(i=0;i<I;i++)
	*(c+i)=*(b+i)/a;
}


//Modified:
void multi(matrix_ &cc,matrix_ &bb,COMPLEX &a)
{
   COMPLEX *c=cc.array,*b=bb.array;
	int i,I=cc.rows*cc.columns;
	for(i=0;i<I;i++)
	*(c+i)=*(b+i)*a;
}
void multi(matrix_ &cc,matrix_ &bb,Complex &a)
{
   COMPLEX *c=cc.array,*b=bb.array;
	int i,I=cc.rows*cc.columns;
	for(i=0;i<I;i++)
	*(c+i)=*(b+i)*a;
}
//solve Ax=b, temp is temporary matrix.
void new_algorithm(matrix_ &A,matrix_ &b,matrix_ &x,matrix_ &temp)
{
	float t3;
   COMPLEX t1,t2;
   int i,j,I,J;
   I=A.rows;J=A.columns;
   for(i=0;i<I;i++)
   {
      t2=COMPLEX(0);
      t3=0;
   	for(j=0;j<J;j++)
      {
      	t2+=b[j]*A[i*J+j];
         t3+=abs(conj(A[i*J+j])*A[i*J+j]);
      }
      t1=t2/t3;
      for(j=0;j<J;j++)
      {
      	temp[j]=t1*A[i*J+j];
      }
      add(x,x,temp);
   }
}
//for two vectors of n degrees
//aH*b;
COMPLEX dot(matrix_ &a,matrix_ &b)
{
	int i,I;
   COMPLEX d=COMPLEX(0);
   I=a.rows*a.columns;

   for(i=0;i<I;i++)
   d=d+conj(a[i])*b[i];
   return d;
}
//<a,b>=bH*a
COMPLEX product(matrix_ &a,matrix_ &b)
{
	int i,I;
   COMPLEX d=COMPLEX(0);
   I=a.rows*a.columns;

   for(i=0;i<I;i++)
   d=d+conj(b[i])*a[i];
   return d;
}
//a*b.
COMPLEX dot0(matrix_ &a,matrix_ &b)
{
	int i,I;
   COMPLEX d=COMPLEX(0);
   I=a.rows*a.columns;

   for(i=0;i<I;i++)
   d=d+a[i]*b[i];
   return d;
}


//This is this the modified H. All the array operations are replaced by pointer operations.
//a is a m-by-n matrix,while b is a n-by-m matrix.   b=H(a)
void H(matrix_ &aa,matrix_ &bb)
{
   int i,j,jn;
   COMPLEX *a=aa.array,*b=bb.array;
   int m=aa.rows,n=aa.columns;
	for(i=0;i<n;i++)
   {
      jn=0;
   	for(j=0;j<m;j++)
	   {
   	   *b=conj(*(a+jn));
      	b++;
         jn+=n;
      }
      a++;
   }
}



//This is the modified duplicate: array operations->pointer operations
void duplicate(matrix_ &aa,matrix_ &bb)
{
  	int i,j;
   COMPLEX *a=aa.array,*b=bb.array;
   int ar=aa.rows;
   int ac=aa.columns;
   int bc=bb.columns;
   for(i=0;i<ar;i++)
   {
   	for(j=0;j<ac;j++)
   	{
   		*a=*(b+j);
      	a++;
      }
      b+=bc;
   }
}

void Tduplicate(matrix_ &a,matrix_ &b)
{
  	int i,j;
   int ar=a.rows;
   int ac=a.columns;
   int bc=b.columns;
   for(i=0;i<ar;i++)
   for(j=0;j<ac;j++)
   {
   	a[i*ac+j]=b[j*bc+i];
   }
}


void duplicate_the_sampling_rows(matrix_ &A,matrix_ &B,int *rowNo,int num_sam)
{
	int i,j,k,l;
   int col=B.columns,r=A.columns;
   k=0;
   for(i=0;i<num_sam;i++)
   {
   	l=rowNo[i]*col;
   	for(j=0;j<r;j++)
   	{
   		A[k]=B[l+j];
      	k++;
      }
   }
}

void imagine_part(float *a,matrix_ b)
{
   int n=b.rows*b.columns;
   for(int i=0;i<n;i++)
   a[i]=imag(b[i]);
}
void real_part(float *a,matrix_ b)
{
   int n=b.rows*b.columns;
   for(int i=0;i<n;i++)
   a[i]=real(b[i]);
}
void foutput(char *f1,char *f2,matrix_ &a)
{
   ofstream fo1,fo2;
   fo1.open(f1);
   fo2.open(f2);
   int i,j,k;
   k=0;
   for(i=0;i<a.rows;i++)
   {
      fo1<<"\n";
      fo2<<"\n";
   	for(j=0;j<a.columns;j++)
   	{
      	fo1<<real(a[k])<<" ";
         fo2<<imag(a[k])<<" ";
         k++;
      }
   }
}
//Compression is used to compress a matrix into a compressed form.
//Thus the new rows and cols should be less or equal to old rows and columns.
void compression(matrix_ &a,int m,int n)
{
   int i,j,k,l,n0=a.columns,m_1=m-1;
   k=n;l=n0;
   for(i=0;i<m_1;i++)
   {
      for(j=0;j<n;j++)
      {
         a[k]=a[l+j];
         k++;
      }
      l+=n0;
   }
   a.rows=m;
   a.columns=n;
}
void Tsquare(COMPLEX &f,matrix_ &a)
{
	int i,I=a.rows;
   f=COMPLEX(0);
   for(i=0;i<I;i++)
   {
   	f=f+a[i]*a[i];
   }
}
void Tmulti(COMPLEX &f,matrix_ &a,matrix_ &b)
{
   int i,I=a.rows;
   f=COMPLEX(0);
   for(i=0;i<I;i++)
   {
   	f=f+a[i]*b[i];
   }
}
void addmulti(matrix_ &aa,matrix_ &bb,COMPLEX c)
{
	int i,I=aa.rows*aa.columns;
   COMPLEX *a=aa.array,*b=bb.array;

   for(i=0;i<I;i++)
   {
   	*a=*a+*b*c;
      a++;b++;
   }

}
void multiadd(matrix_ &aa,COMPLEX c, matrix_ &bb)
{
	int i,I=aa.rows*aa.columns;
   COMPLEX *a=aa.array,*b=bb.array;
   for(i=0;i<I;i++)
   {
   	*a=*a*c+*b;
      a++;b++;
   }

}
void submulti(matrix_ &aa,matrix_ &bb,COMPLEX c)
{
	int i,I=aa.rows*aa.columns;
   COMPLEX *a=aa.array,*b=bb.array;
   for(i=0;i<I;i++)
   {
   	*a=*a-*b*c;
      a++;b++;
   }

}      */

#endif